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Parametrically modulated optomechanical systems have been recently proposed as a simple and efficient set- 
ting for the quantum control of a micromechanical oscillator: relevant possibilities include the generation of 
squeezing in the oscillator position (or momentum) and the enhancement of entanglement between mechanical 
and radiation modes. In this paper we further investigate this new modulation regime, considering an optome- 
chanical system with one or more parameters being modulated over time. We first apply a sinusoidal modulation 
of the mechanical frequency and characterize the optimal regime in which the visibility of purely quantum ef- 
fects is maximal. We then introduce a second modulation on the input laser intensity and analyze the interplay 
between the two. We find that an interference pattern shows up, so that different choices of the relative phase be- 
tween the two modulations can either enhance or cancel the desired quantum effects, opening new possibilities 
for optimal quantum control strategies. 



I. INTRODUCTION 

Theoretical studies and huge technological progresses over 
the last decades made it possible to reach a considerable level 
of control over quantum states of matter in a large variety of 
physical systems, ranging from photons, electrons and atoms 
to bigger solid state systems such as quantum dots and su- 
perconducting circuits. This opened the possibility for novel 
tests of quantum mechanics and allowed, among other things, 
to take important steps forward in investigating the quantum 
regime of macroscopic objects. In this perspective, one of the 
main goals in today quantum science is controlling nano- and 
micromechanical oscillators at the quantum level. 

Quantum optomechanics |'T'-'4'|, i.e. studying and engineer- 
ing the radiation pressure interaction of light with mechani- 
cal systems, comes as a powerful and well-developed tool to 
do so. First, radiation pressure interaction can be exploited 
to cool a (nano)micromechanical oscillator to its motional 
ground-state Q; this is a necessary step for quantum manip- 
ulation and could not be accomplished by direct means such 
as cryogenic cooling (at the typical mechanical frequencies 
involved IQOKHz ^ IGHz this would require cooling the 
environment to a temperature of the order l/iJ-iT ^ lOraK). 
Backaction cooling has been experimentally demonstrated for 
a variety of physical implementations, including micromirrors 
in Fabry-Perot cavities [6J, microtoroidal cavities [7J or op- 
tomechanical crystals O. Second, there exists a strong anal- 
ogy between quantum optomechanics and non-linear quan- 
tum optics, so that many (if not all) optomechanical effects 
can be mapped onto well-known optical effects. As a result, 
optomechanics becomes a natural way for controlling a me- 
chanical resonator at the quantum level. Experimentally, the 
strong coupling regime needed to observe quantum behaviors 
has been demonstrated only very recently [7^^, and detection 
of quantum effects is still awaiting. Nevertheless, a lot of the- 
oretical studies on the subject has been carried out in the last 
decade and several proposals have been produced [ 10|. These 
cover, among other things, the generation of entanglement be- 
tween one oscillator and the radiation in a Fabry-Perot cavity 
ifTTl . the generation of entanglement between two oscillators 
lfT2l . or the generation of squeezed mechanical states il3l[T4l . 
In particular, references lfT4HT6l introduced a new and ef- 



fective way of enhancing the generation of quantum effects, 
which relies on applying a periodic modulation to some of the 
system parameters (a similar result has also been found in the 
analogous contest of nanore senators and microwave cavities 
|17J). 

In this paper we further investigate the properties of peri- 
odically modulated optomechanical systems and we address 
the following questions: which is the fundamental link be- 
tween modulation and enhancement of quantum effects? is 
there an optimal choice of the modulation, for which the vis- 
ibility of quantum effects is maximal? Is this optimal regime 
robust against parameter fluctuations? What happens when 
two independent modulations are applied simultaneously? To 
tackle these issues we analyze the paradigmatic case of a me- 
chanical oscillator whose natural frequency lom is externally 
modulated when it evolves under the action of the noise and of 
the radiation pressure exerted by the photons of an externally 
driven optical cavity mode. While quantum optomechanics is 
nowadays extensively studied within a variety of experimental 
setups, the modulation of the mechanical frequency we ana- 
lyze here is a very crucial aspect of our system and one that 
has not been implemented yet. However, very recent propos- 
als for doing optomechanics with levitated dielectric spheres 
ifTSl - EOl can be a good answer In these proposals the mechan- 
ical degree of freedom is represented by the center of mass 
motion of a nanodielectric sphere which is trapped and lev- 
itated by means of an optical trap. The sphere is then put 
inside an optical cavity, where it interacts with the intracavity 
radiation via the usual optomechanical Hamiltonian ([2]). The 
frequency of the center of mass motion depends on the shape 
of the trapping potential and can thus be modulated adjusting 
the intensity of the trapping laser, as shown in 1 1 8 1 . Moreover, 
typical parameters attainable with such setups are comparable 
to those we have adopted in our simulations (see below), as- 
suring the feasibility of the system under analysis. 

In the above scenario we study the formation of squeez- 
ing, entanglement and discord f2r|, showing that in the steady 
state all these quantum effects are enhanced when the mod- 
ulation frequency Q. is twice the original value of lom- As 
we shall see, such resonance admits a simple inteipretation in 
terms of an effective parametric phase-locking between the 
external driving forces and the natural evolution of the in- 
volved degree of freedom. Similar enhancements were also 
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Figure 1 . Schematic description of the system. A Fabry-Perot cav- 
ity is driven by an external laser and the radiation interacts with the 
movable mirror on the right, exchanging momentum. 



observed in Refs. lfT4l[T5l . where an harmonic modulation of 
2lom was imposed on the amplitude of the cavity mode laser, 
and in Ref. |16 |, where a harmonic modulation of 2liJm was 
imposed on the coupling rate between two generic bosonic 
modes. Since several mechanisms can lead independently to 
the same effect, an interesting question is how they can be best 
exploited to control specific quantum properties in the system. 
This goes in the direction of developing optimal quantum con- 
trol protocols, a topic which is currently benefiting from many 
contributions |22|. In the present case, to study the interplay 
of different mechanisms we add a second modulation in our 
model and we observe the arising of interference pattern in 
the system response. Specifically we notice that the ability 
in cooling and squeezing the mechanical oscillator strongly 
depends upon the relative phase of the two modulations, the 
relative variation being almost 50%. 

The material is organized as follows. In section |ll] we 
present the system and solve its dynamical evolution under the 
action of a periodic modulation of the mechanical frequency. 



In section III we then characterize the asymptotic stationary 
state in terms of entanglement, squeezing, etc. In sectio n |IV| 
we compare our findings to other recent proposals fl4l 1151 
and we study what happens when a second independent mod- 
ulation is applied to the system [specifically, in our case we 
introduce a modulation on the amplitude of the input laser]. 
Conclusions and general remarks follows in section [V] S ome 
technical derivations are finally reported in Appendix |A| 



II. THE SYSTEM 

Our choice falls on the simplest optomechanical system of 
all, i.e. a Fabry-Perot cavity of lenght /q with a movable mirror 
at one end (see Fig. [T]i, which nevertheless captures all inter- 
esting physics. We can reasonably assume |10| that a single 
optical mode is interacting with a single mechanical mode, 
be it the center of mass oscillation. The mirror can thus be 
modeled as a mass m attached to a spring of characteristic 
frequency uj^ and friction coefficient 7m; it is described by 
dimensionless position and momentum operators q, p which 
obey the canonical commutation relation [q,p\ = i. The op- 
tical mode has frequency uc and decay rate k; it is described 
by annihilation/creation operators a and a\ which obey the 
canonical commutation relation \a,a)] = 1. 



In our analysis the cavity is assumed to be driven by an 
external laser which, to begin with, we take to have constant 
power Piaser ^nd quasi-resonant frequency oj^ ^ ljc- In this 
context a periodic modulation is inserted at the level of the 
spring constant, which we express as the following time de- 
pendent parametric rescaling of the mirror frequency 



[1 + ecos(m)], 



(1) 



with e < 1. Accordingly the Hamiltonian of the system writes 
as Gl 



hGoo} aq + ihE{e 



(2) 



where Go = i^c/h i'm^M) is the optomechanical cou- 
pling rate and \E\ = -^Z 2 fc Piaser /^l is the driving rate. In- 
cluding dissipation and decoherence effects the system dy- 
namics can then be described with the following set of quan- 
tum Langevin equations ifTOl 

dt q = iOMP, 

dtp= -ujm [1 + e cos(rit)] q - jmP + Gao^a + ^, 
dt a = — (fc + zAo)a + iGoaq + E + y/2k aim 

(3) 

which we have written in a frame rotating at ujl. Here 
Aq = ujc — is the unperturbed cavity laser detuning while 
ain{t) is the radiation vacuum input noise with autocorrela- 
tion function I 



{t)al{t'))^5{t~t'). 



(4) 



Similarly ^(t) is the Brownian noise operator describing the 
dissipative friction forces acting on the mirror. Its autocorre- 
lation function satisfies the relation ||25| 



(5) 

which for the specific case of an harmonic oscillator with a 
good quality factor ojm ^ 7Af, acquires the same Markov 
character of Eq. Q, i.e. 

[mA{t')]) « 27MCoth(^) 8{t-t'), (6) 

(this is a consequence of the fact that for lom ^ jm only 
resonant noise components at frequency u! ^ lum do sensi- 
bly affect the motion of the system). In the above expres- 
sions T is the system temperature while {• • • , • • • } is the anti- 
comutator L26j. 



A. Solving the dynamics 

The evolution of the system is ruled by a set (|3]l of non- 
linear stochastic differential equations with periodic coeffi- 
cients, whose solution is in general very difficult. In the fol- 
lowing we will then introduce some useful approximations to 
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simplify the calculations. First, we expand each operator as 
the sum of a c-number mean value and a fluctuation operator, 
i.e. 



a{t) = (a(t)) + {a{t) 



(7) 



We recall that the cavity is usually driven by a very strong 
laser in order to attain satisfactory levels of optomechanical 
interaction, so that the mean value will be much bigger than 



the fluctuations, which are due to the presence of random 
noise. This allows us to write (p) as two different sets of equa- 
tions, one for the mean values^S}, one for the fluctuations ^ 
and linearize the latter neglecting all terms which are second 
order small, obtaining 



dtQ = ujAiP, 

dtP^ -UJM [1 + ecos(m)] Q - -imP + Go\A\\ 
dtA^ -(fc + iAo)^ + iGaAQ + E, 



(8) 
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(9) 



where we have introduced the phase and amplitude quadra- 
tures for the cavity and the input noise fields, i.e. 

X = jat+a)/y2, Y - i{a^-a)l^, X,„ - (aL+«m)/^ 
and Yin ~ ^('^L ~ am)/\/2- Equation can be also ex- 
pressed in a more compact form 



dt u — Sit + C, 



(10) 



with S being a 4 x 4 time-dependent matrix, and with u and 
( being the column vectors of elements {Sq, 6p, SX.SY) and 
(0, 1, Xin, Yin), respectively. We stress that Eqs. (jsjl and ^ 
must be solved in the correct order, because the mean values 
Q{t), P{t) and A{t) play the role of coefficients in the equa- 
tions for the fluctuations. 

Equation ([8]l is nonlinear but can be solved numerically. 
Assuming that we are far from optomechanical instabilities 
and that we keep the modulation strength e small enough to 
avoid additional instabilities due to parametric amplification, 
one finds that the mean values evolve toward an asymptotic 
periodic orbit with the same periodicity 27r/ri of the applied 
modulation. In this regime, an approximate analytic solution 
can also be derived, which we detail in Appendix |A] Indeed 
since the modulation strength e is not too strong, one can guess 
a perturbative expansion of the form 



g(t) = ^Q(^")(i), 



(11) 



and the noises are zero-mean gaussian noises, fluctuations in 
the stable regime will also evolve to an asymptotic zero-mean 
gaussian state. The state of the system is then completely de- 
scribed by the correlation matrix G of elements 

G^j{t) = Gj,{t) = i {U;{t)Uj{t) + Uj{t)u^{t)) , (12) 



whose evolution can be derived directly from equations ( lOl 
and ([T2]i: 



dtG = SG + GS^ 



(13) 



where is transpose of S, and where N is the 

diagonal noise correlation matrix with diagonal entries 

(0, 7m coih{hjjM /'i-ksT), fc, k), defined by 



C.(t)0(t') + Q{i')Ut)) = N,j5{t - t'). (14) 



Equation ( 13 i is now an ordinary linear differential equation. 
We know that its solution evolves toward a unique asymptotic 
configuration (independently of the initial state), proven that 
the eigenvalues of the matrix S have negative real part for all 
times t, which can be verified by applying the Routh-Hurwitz 
criterion lIZTll . Again, we can either solve Eq. (13 1 numerically 
or obtain an approximate analytic solution with a perturbative 
expansion in e (see Appendix [A] for the latter). 



where Q^^^t) does not depend on e, Q'^^^t) is linear in e, 
Q^'^\t) is quadratic in e and so on. It turns out that each or- 
der is exactly solvable, as long as previous orders are known. 
This originates a chained set of equations and by keeping a 
finite number of orders j < Jmax, we can finally obtain the 
asymptotic solution up to the desired precision (e.g. see Fig. 

Equation (j9]l is stochastic and needs some more manipu- 
lation. Nonetheless since we have linearized the dynamics 



B. Quantum properties of the system 

As already mentioned, thanks to gaussianity of the asymp- 
totic solution all relevant informations about the system can 
be extracted directly from the correlation matrix G. In par- 
ticular we will focus on the following quantities: the number 
of phonons in the mirror, the squeezing in the mirror and in 
the radiation quadratures, and the nonclassical correlation be- 
tween the mirror and the radiation degrees of freedom. 
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Analogous formulas for the radiation quadratures lead to 
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Figure 2. Evolution of the mirror position Q{t) and momentum mean 
values P{t) obtained by numerically integrating Eq. (jsjl from t — 
to t — 50r, with r = 2ti/Q. being the period of the modulation 
(thin blue line). The plot has been obtained by setting the system 
parameters as detailed in Sec.|Ill] in particular here the modulation 
frequency Q, is twice the natural frequency ujm of the mechanical os- 
cillator which, in turn, is resonant with the detuning Ao that governs 
the free evolution of the optical field A{t). The analytic solution for 
the asymptotic orbit (see Appendix ) A 3) is also shown for compari- 
son (thick red line). 



The number of phonons n can be expressed using the ap- 
proximate relation 



hujul^ (Cii 



C12) , 



which holds if the modulation of the mechanical frequency 
is not too strong. This tells how far the system is from the 
ground state. Since both Cn and C12 are periodic in time, we 
will identify the number of phonons with the maximum over 
one period r = 27r/f7 of the modulation, i.e. 



nuAX 



:Mt)}, 



(16) 



(here and in the following max^ represents an optimization 
with respect to a time interval [T, T + r] with T being suf- 
ficiently larger than 1/fc to guarantee that the system has 
reached the asymptotic steady state). 

Squeezing of the generalized mirror quadratures qs — 
q cos 9 + p sin 9 is also easily found: 

(Sqj) = Cii cos^ 9 + C22 sin^ 9 + (Ci, + C21) cos 9 sin 61. 

(17) 

Again we construct a time independent quantity to deal with. 
First, for each time t we select the parameter 9 for which 
{Sqg') is minimum. In terms of covariance matrix, this is just 

the smaller eigenvalue of the block matrix ( „ „ 

then minimize this quantity with respect to time over a period 
T. This tells how much squeezing can be produced at most. 



A'^QAfjN = min <^ min (dqg) 



(18) 



MIN 



mm < mm 



(19) 



Non-classical correlations in the system can be described us- 
ing quantum discord [21], which includes entanglement 
as well as more general quantum correlations that are shown 
also by separable states f2E\. For a gaussian state, is eas- 
ily constructed from the correlation matrix as demonstrated in 
II29II . Time dependance is then eliminated by considering 



'Dm AX = max{X'(/9(t))} , 



(20) 



Entanglement alone will be specifically described using log- 
arithmic negativity Ej^{p) lf30l . which is also easily con- 
structed from the correlation matrix as demonstrated in |31|. 
Again, time dependance is eliminated by considering 

E^MAx ^^i^^iEMpim- (21) 



III. RESULTS 

We now present the results obtained by solving the 
dynamics of the system as detailed in the previous section. 
The parameters used in our analysis are m — 150 ng, 
u}m/{2t:) = 1 MHz, -fM/{2n) = 1 Hz, T = 0.1 K, 
Ao = LJM, lo = 25 mm. A: = 1.34 MHz, A = 1064 nm and 
flaser = 10 mW: this choice is compatible with values attained 
in state of the art experiments and is also consistent with the 



(15) stability requirement of section II A (furthermore, under the 



condition Ao = ujm the optical and the mechanical variables 
are brought at resonance). The strength e and the frequency 
51 of the modulation are left as variable parameters instead, 
since we want to characterize the optimal modulation regime, 
e.g. which e and fl maximize the visibility of quantum effects. 

In Fig. |2] we temporarily fix 51 = "^com, e ~ 0.2 (this 
particular choice will be justified in the following) and we 
report the solution of Eq ^ for the mean values Q{t) and 
P{t) of the mirror position and momentum. We see that the 
evolution tends indeed to an asymptotic periodic orbit, which 
is very well approximated by the analytic solution. 

We then focus on the solution of equation ( [T3| ) and we plot 
the quantities described in section II B for multiple values of 
fl G [(jJm e G [0,0.5]. In particular: Fig. [sjshows the 
maximum number n max of phonons in the mirror, computed 
via Eq. ( [T6] i; the maximum ^JaAmay °f '^^e logarithmic nega- 



tivity, computed via Eq. (21 1; the maximum Vmax of the 



quantum discord, computed via Eq. (20i; and the minimum 
variance A'^qMiN of all the mirror generalized quadratures, 



We computed via Eq. ( 18 i 



As evident from the plots, the level of squeezing and en- 
tanglement is maximum when the modulation frequency is 
51 2u!m and increases monotonically with respect to the 
strength e, until the system eventually reaches an instability 
point for too strong modulations (in the above figures, this 
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instability is represented by a blank region around the point 
e = 0.5, il ~ 2ujm)- It is also clear that the optimal mod- 
ulation, the one that most enhances quantum effects, is also 
responsible for heating the system far from its ground state. 
We can understand this behavior if we interpret Eq. (13i as 



describing the dynamics of a set of (classical) parametric os- 
cillators with canonical coordinates defined by the correla- 
tions functions (12i, which evolve under the action of 
damping and constant external driving forces. Indeed, by a 
close inspection of the matrix S one notices that such os- 
cillators possess natural frequencies which are periodically 
modulated through functions (i.e. A{t), Q{t), and the di- 
rect term ecos(i7t)) that, in first approximation, evolve sinu- 
soidally with the same frequency fl - see Eq. (A17i in Ap- 
pendix |A] for details. Moreover, in the stability region we are 
sure that parametric modulation pumps energy into the sys- 
tem at a lower rate with respect to losses, since the system 
evolves toward a stationary orbit: we call this regime "below- 
threshold" to distinguish it from the exponential amplification 
usually associated with parametric oscillators. For this model 
phase-locking is expected to occur when fl matches the zero- 
order eigenfrequencies defined by the constant part of S (and 
not twice this frequencies as in the case of parametric instabil- 
ity), resulting in an enhancement of the oscillations of the ef- 
fective coordinates Cij ( [T2| and hence of the associated quan- 
tum effects defined in Sec. |IIB| [321 (more details are found in 
Appendix [A|. It turns out that, at least for the figure of merit 
we are concerned here (i.e. A'^q]\fjN, I^^Xmin, T^max, etc) 
the relevant frequency is indeed ^ Iujm- 

To see this, we can procede by steps. First of all notice 
that from the numerical solution, we can guarantee that the 
system is not unstable (see Fig. [3]l, i.e. that it is indeed in 
the below-threshold regime. Next, consider the case of no 
coupling (Go = 0) and no modulation (e = 0): we stress 
out three relevant aspects. First, the mechanical part and the 
radiation part are independent, so there is no entanglement. 
Second, each subsystem evolves with the Hamiltonian of a 
quantum harmonic oscillator, so the quadrature mean value 
{qe) evolves with a phase e*'^"* and the variance ((J^g) with 
a phase e'^"^*'* (we remind that we fixed Ag = wa/)- This 
tells us that, at least in this regime, the frequencies which gov- 
ern the quantities of interest are degenerate at the value 2ljJm- 
Third, each subsystem is also coupled to its own environment 
and will eventually relax to a thermal state characterized by 
{^Qs) ^ ^ therm + 1 /2, SO there is no squeezing. Now turn on 
the coupling Gq: this has three main effects. First it introduces 
entanglement in the system 1,1 1 J {Ej^/ — Eq). Second the 
eigenfrequencies are brought out of degeneracy and shifted 
by a term cx 2Go \ A\ (33], which is quite small with respect to 
2ujm for our choice of values (confirming that indeed the latter 
is the resonant value at which the modulation should provide 
an enhancement). Third, backaction cooling |5 | is now active 
and the oscillator approaches the ground state {{Sqg) ^ 1/2). 
Squeezing is still absent at this level. Finally, turn on the mod- 
ulation ([T]). Thanks to the phase-locking mechanism we have 
anticipated previously and detailed in Appendix [A] this will 
yield an enhancement of the correlations when 51 matches the 
natural frequency ^ 2ujm- For instance, for the negative en- 



tropy Ef^ and for mirror variance (Sqg), we get 

E^f ^ Eo + e Ki{n) cos{nt + (pi), 
(H')o ^ 1/2 + e K2{n) cos{nt + ^2), 



(22) 



where Ki{n) and K2{i^) are associated response functions 
analogous to the Lorentzian response of a simple harmonic 
oscillator (though an exact expression is rather cumbersome 
in our specific case) and are peaked around ft ^ 2ujm- We 
see that the quadrature Sq^ gets periodically squeezed over 
time and entanglement is periodically increased to higher val- 
ues with respect to the unmodulated case. In addition, these 
effects increase monotonically with e, up to the instability 
threshold. A similar enhancement of the entanglement is also 
described in Ref. III6I . where two harmonic oscillators are 
coupled via linear interaction Hint = cit)XiX2 and the cou- 
pling constant is a periodic function of time. This time de- 
pendance produces an effective modulation on the normal fre- 
quencies of the system: as a result, entanglement is shown to 
increase and become much more robust against temperature. 
This agrees very well with what we found here. 



IV. INTERPLAY BETWEEN TWO DIFFERENT 
MODULATIONS 

Results analogous to those presented in the previous sec- 
tion has been found very recently by Mari and Eisert [14, 15|, 
for an optomechanical system driven with an amplitude mod- 
ulated input laser. For clarity, we rewrite their Hamiltonian 

"H =fujjca}a H ^ (p^ + g^) — hGoo'aq 



ih{E + Ei cos(m))(( 



a 



(23) 



At first sight, the situation appears to be somewhat different 
from our initial problem. In Eq. ( 23 1, internal parameters of 



the system are left unchanged; it is instead the external driv- 
ing that undergoes an oscillatory behavior. Nevertheless the 
effects are strikingly similar: high levels of squeezing can 
be attained when the frequency of modulation is ^ "^ujm 
IIT4II . and the same regime is also optimal to enhance entan- 
glement between mechanical and radiation modes |15| The 
authors themselves comment that "...this dynamics reminds 
of the effect of parametric amplification, as if the spring con- 
stant of the mechanical motion was varied in time with just 
twice the frequency of the mechanical motion, leading to the 
squeezing of the mechanical mode..." lfT4l l. 

In fact there is a strong analogy between the two cases. In- 
dependently of what Hamiltonian (j2|i or ( [23] l one chooses, far 
from instability regions the mean values Q{t), P{t), A{t) will 
be characterized by an asymptotic periodic orbit with the same 
periodicity of the applied modulation t = 27r/il. This assures 



that in both cases the equation ( 13 1 for the covariance matrix 



has the same linear form, with S being a periodic function of 
time (in the limit i ^ 1/fc) and N being a constant driving. 
The conclusions we derived in section|lIl] must therefore hold, 
at least qualitatively, also for the system studied in llT4irT5]| . 
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Figure 3. Asymptotic quantum features as a function of fl/uiM (x axis) and e (y axis): a) Maximum number of phonons in the mirror, 
eq \16\ ; b) Maximum of the logarithmic negativity \2l\ ; c) Maximum of the quantum discord |20|; d) Minimum of the generalized 
quadratures of the mirror JTSl. In all the plots the system parameters are fixed as in Sec.|III| 



An interesting question now rises. What if the two modula- 
tions are applied together? Can they interfere, either construc- 
tively or destructively, and sensibly alter the one-modulation 
picture? 

To get an answer, we consider a new composite system, 
described by the Hamiltonian 



H = huca' a H —p 



2 



1 



ecos 



-hGoo'aq (24) 
+ihE{l + ri cos{n2t + 0))(e-*'^^*a^ - e''"^*a). 



Note that we explicitly introduced a relative phase between 
the two applied modulation; if we expect any interference, the 
properties of the system should indeed depend on this new 
variable. 

The analysis presented in the previous sections is straight- 
forwardly generalized to the present case, so we will skip di- 
rectly to the results [details can be found however in the Ap- 
pendix 



III 



Taking the same parameters as in Sec 
the optimal modulation frequencies fli = fl^ = 
fix e = 0.3, f] = 0.9 (this is the same value used in 



we choose 
2ujm and 
141). 



These modulation strengths give comparable squeezing per- 
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Figure 4. Response of the system in the presence of two different modulations as a function of their relative phase 0: a) Maximum 
number of phonons in the mirror; b) Maximum of the logarithmic negativity; c) Maximum of the quantum discord; d) Minimum of the 
generalized quadratures of the mirror. In all the plots the two straight lines show the variation of the function in the case when only the 
mechanical frequency (blue horizontal) or the laser amplitude (red dashed horizontal) is modulated. Parameters as detailed in the text. 



formances when considered singularly, and also assure that we 
are reasonably far from the instability region. To present the 



results, we plot the quantities introduced in Sec. II B against 
the relative phase </) in Fig. |4] An interference pattern is in- 
deed evident and each of the above quantities oscillates be- 
tween a minimum and a maximum as cf) varies in the range 
[0, 2tt]. However, entanglement and quantum discord are af- 
fected very weakly and do not differ much from our initial 
one-modulation case. Besides we see that in order to gen- 
erate quantum correlations, a modulation of the mechanical 
frequency is more suitable than a modulation of the driving 
laser amplitude. Adding the second modulation to the first is 
of little effect. 

Squeezing generation instead, presents very interesting fea- 
tures. First, as we said, we choose two modulations that 
give comparable levels of squeezing when applied individ- 
ually. Moreover, when applied together, they can strongly 
interfere. For example we see in Fig. [4|d) that for a phase 
(j)/TT ~ 1.4, A'^qmin rises toward the threshold value 0.5 
and squeezing becomes weaker. Each modulation taken alone 
would generate more squeezing than the two combined: this 
is an unambiguous sign of a disadvantageous interplay. For a 
phase (j>/Tr ^ 0.4 we find instead a great advantage in apply- 
ing two modulations: A'^qmin is lowered to a value ^ 0.18, 
a considerable performance if compared to our initial one- 
modulation case where instabilities prevent us from reaching 
A'^qmin < 0.17. In fact, not only we attain the same high 



levels of squeezing, but we are also well inside the stability 
region, so that we could increase both e and rj to perform even 
better. 

We also note that the optimal(worst) phase choice for 
squeezing generation also corresponds to maximum heat- 
ing(cooling) of the mirror, as can be seen in Fig. |4|^)- This is 
another confirmation that parametric oscillation is indeed the 
main underlying mechanism: in fact not only does a stronger 
modulation enhance the generation of quantum effects, as in- 



ferred from equations ( 22 1, it also pumps more energy into the 
system. 

We then see how the interplay between two independent 
modulations can be carefully exploited to increase levels of 
squeezing in an optomechanical system. 



V. CONCLUSIONS 

We have studied in great detail the effect of periodic mod- 
ulations on optomechanical systems and we have character- 
ized several ways in which such modulations can be exploited 
to enhance relevant quantum properties including squeezing, 
entanglement and quantum discord. While the idea that mod- 
ulations can help accessing the quantum regime was already 
known from previous works 1 14- 16 1, we have proposed a new 
interpretation of this enhancement mechanism in terms of a 
resonance between the modulation frequency and the natural 
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frequencies of the system. This simple model allowed us to 
prove the existence of an optimal modulation regime and to 
understand the arising of instability thresholds. Finally, we 
have analyzed the interplay of different modulations and we 
have found that constructive (destructive) interference effects 
may arise when they are applied simultaneously, causing a 
further enhancement (a suppression) of quantum effects. We 
beheve that these results could lead further on toward the de- 



velopment of optimal control strategies. 
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Appendix A: Asymptotic behaviors 

This section deals with some technical aspects related with 
the asymptotic solutions of Eqs. (l8]l and ( 13 1, which define the 



quantum properties of the system. Here we discuss the reso- 
nant mechanisms which is responsible for the enhancement 
of quantum effects at 57 ^ 2cl>m, as well as the role of the 
relative phase in the interplay between different modulations. 



We start in Sec. A 1 by presenting a simple paradigmatic case 



which captures the main aspects of the resonance. Then, in 
section |A 2| we introduce the analytic framework which will 
be used to describe the dynamics of the system. Finally, sec- 
tions A 3 and A 4 are devoted to analyze in details the asymp- 



totic behavior of the system, in the one- and two-modulation 
scenario respectively. 



1. Single oscillator model 



with f{v) = (wq — v^)^ + {ji^)^ being the Lorentzian 
response function of a classical harmonic oscillator Clearly 
the superimposed oscillation, which we remind is an effect of 
the parametric modulation, will be much greater near reso- 
nance with the natural frequency ly ^ ojq and for a just below 
the instability threshold. Going to second order in a yields 
small deviation from this picture and we can stop our quali- 
tative analysis here. In summary, parametric modulation can 
controllably enhance oscillations of the system coordinates if 
two main conditions are satisfied: the modulation must not 
be too strong, otherwise the system becomes unstable, and 
an external (constant) driving must also be applied, other- 
wise the system relaxes to x{t — >■ oo) = (as from Eq. A3 
with F = 0). We also stress out that, in the below-threshold 
regime, the resonance condition is given hy ^ loq (i.e. the 
modulation frequency should be the same as the natural fre- 
quency of the system) and not by '-^ 2ujo, as is the usual case 
of exponential parametric amplification. 



As anticipated in the main text the evolution of the corre- 
lations matrix describes the dynamics of a multi-dimensional 
(classical) oscillator which evolves in presence of damping 
and external constant driving (defined by the matrix N) and 
which possesses characteristic frequencies (determined by 
S{t)) that are externally modulated at frequency 51 [these 
statements are explicitly verified in Sees. |A 2] |A 3| and |A 4| . 
To enlighten the role of the modulation in the evolution of the 
correlation functions it is hence worth focusing on the sim- 
plest example of this sort. This is provided by a single para- 
metric oscillator whose position x evolves according to the 
equation 



x{t) = -cjo[l + acos{vt)] x{t) - 7i(i) + F 



(Al) 



with uiQ and being the characteristic and the modulation fre- 
quency, a being the amplitude of the modulation, 7 being the 
damping rate and F the strength of a constant driving. For this 
simple scenario, two cases are possible. If a > (above- 
threshold condition), parametric modulation pumps energy 
into the system at a faster rate with respect to dissipation; the 
system increases its energy exponentially and is therefore un- 
stable. If a < 27/^0 (below-threshold condition), the system 
reaches a stationary regime, given by the balance of pump- 
ing and dissipation. We can then look for a stable solution 
of Eq. ( Al I, assuming that a is small and treating the solu- 
tion perturbatively, i.e. x{t) = + ax^-'-^t) + 0{a^). 
To order zero in a the system is just a damped driven har- 
monic oscillator, which relaxes toward its equilibrium posi- 
tion a;^°) = limt_j.oo x^^'>{t) — F/oj'^. To first order in a, the 
long time solution is then given by 

-ujlx'^^\t)~u^cos{iyt) -7i(i)(i). (A2) 

Therefore we see that the parametric modulation, for the 
below-threshold regime, can be mapped onto an effective ex- 
ternal driving F cos(i^t) and the solution is easily found to be 



F 



2. General treatment of the modulated optomechanical system 

Turning back to Eq. ( [T3] l, we will see that all conditions 
are indeed satisfied: the coefficient matrix S{t) is periodically 
modulated over time, stability can be verified with a numeric 
solution and external driving is provided by the noise corre- 
lation function N. The above result implies that in the case 
of our multi -dimensional parametric oscillator, maximum en- 
hancement of the oscillations is expected when ft matches the 
characteristic frequencies that govern the dynamics of the cor- 
relation functions in absence of the modulation. The latter are 



defined by the matrix S{t) of Eq. (13 1 when e = (and rj = 
in the two-modulation scenario). As mentioned in the text, 
at least for the figure of merit we are concerned about in the 
paper (i.e. A'^qmin, A'^Xmin, T^max, etc), the relevant 
frequency is indeed ^ 2lum- 

We can thus generalize the simple model of Sec. A 1 to the 
present case and reproduce the numerical results we found in 
the main text with a semi-analytic solution of Eqs. (|8]l and 
( 13 1, which we briefly sketch here. In doing so, we will also 
identify and comment the relevant points which are responsi- 
ble for the behaviours observed in Figs. |3]and|4] 



a. Classical solution. 

Let us start with Eq. ([8]l for the mean values which, for the 
sake of completeness, we report here for the general scenario 



a f{i') F cos{i't + I 



(A3) 



defined by the Hamiltonian ( 25 1 where both the frequency 
modulation (j2]) and the amplitude modulation ( 23 1 are acti- 
vated, i.e. 

dtQ ^ ujmP, 

dtP^-UM [l + ecos{nt)]Q-^MP + GM\\ 
dtA = -{k + iAo)A + iGoAQ + E[l + f]cos{nt + 0)], 

(A4) 

having only assumed their frequencies to be identical, i.e. 
ill = ^l2 = ^- As anticipated in the text - see Eq. ( 1 1 1 - we 
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look for a perturbative solution in the modulations strengths e 
and 7], i.e 



Q = QW + 0(i) + g(2) + ..., 

^ = Aio) + + + . . . , 

where, for instance, Q*^^'' is linear in e and 77, Q^^^ is quadratic 
in e and 77 and so on (note that we can revert to the single 
modulation scenario simply by imposing rj = 0). At order 



zero we get 

< dt P(o) - -COM Q(") -7A/P(") +Go|A(o)|2, 

(A6) 

From the numeric simulation we know that this non-linear 
equation evolves toward a stable point (Q^^K P''^\ A^^^^) 
and by setting the derivatives to zero, we can find these 
asymptotic values. Next, at first order we get 



dt P(i) = -UM 0^1) - OJM ecos{nt) Q(o) - 7mP(i) + Go A^i) + Go (^(i))* (A7) 

dt = -{k + iAo)A(i) + iGoi(")Q(i' + iGoA(i)Q(o) + Er] cos{nt + </.). 



These are the equations of three coupled and forced harmonic 
oscillators, with forcing terms ~ujm ecos(ilt) Q'-'^' and 
E rycos(ilt + cf)) that are purely oscillating. In addition, 
damping makes sure that the system is stable. The solutions 
are easily obtained in the form 



with qi, pi, ai, and 02 being complex parameters which can 







= qieJ-^*- - 






1 


p(l)(<) 






(A8) 




AW(t) 


= aie'^^^ - 







be computed by replacing ( A8 1 into Eq. ( A7 1. Hence, to first 
order in e and rj, the effect of the modulation on classical 
values is to add an oscillating term with frequency ft and 
mean value 0. The amplitude of this oscillation clearly de- 
pends on the forcing term, i.e on the amplitudes e, 77, on their 
relative phase (p and on the frequency 51 (via the oscillator 
response function). Finally, the equations for second order are 



'dtQ('^ =lomP^'\ 

dt P(2) = -ijM Q^'^ - ecos(m) Q(i) 
dt = -{k + iAo)A(2) + iGoA(2)Q(o) - 



' 7mP(2) + GoA*(")A(2) + GoA^^^AW + GoA*(2)^(o)^ 
^GoA(l)QW+^Go^(")0(2)_ 



(A9) 



These are again the equations of three coupled and forced 
harmonic oscillators, but this time the forcing terms 
-OJM ecos(m) Q(i), GoA*(i)A(i) and iGoA^^^Q'^^^ have 
also a constant part. As for the first order corrections the 
solutions are easily obtained in the form 

< p(2)(i) = p(2) +p3e^2nt _^p*g-»2nt^ (^10) 

To second order in e and 77, the effect of the modulation on 
classical values is thus to add a constant shift and an additional 
oscillating term with frequency 257 and mean value 0. Higher 
orders can be processed in the same way but for the parameter 



region we have selected in the main text, one can limit the 
analysis to second order since already at this point we get the 
correct result within a good degree of accuracy (see Fig |5]l. 
Full convergence of the approximation when higher orders are 
included can be seen from Fig |2] in the main text, where we 
plot the numeric evolution of classical values Q and P and 
the analytic counterpart, computed up to order six. 

b. Linearized quantum solution. 

We now turn to Eq. (j9]l for the quantum fluctuation, which 
we rewrite below 

dtC = SC + CS^ + N. 
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Figure 5. First order (red dashed) and second order (blue solid) ap- 
proximation to the classical position-momentum (Q-P) orbit of the 
mirror. The asymptotic numerical orbit (black dotted) is also plotted 
from t — 49t to t = 50r, with r = 2-k/Q, being the period of the 
modulation. It is clear that a first order approximation is not enough 
and fails to describe the dynamics of the system. On the other hand, 
the second order approximation catches all relevant aspects and re- 
produces the correct behavior, at least on a qualitative level. Quanti- 
tative convergence to the numeric solution is found including higher 
order terms. 



C = C(")+C(i)+C(2) + .... (All) 

where C^^^ is linear in e and r], (7^^^ is quadratic in e and rj 
and so on. The calculations simply follow what we have done 
for the classical part. At order zero we get 

dt = • + • + N. (A12) 

This equation is linear and evolves toward a stable point C'('^\ 
which we can find by setting the derivatives to zero. Next, at 
first order we get 

(A13) 

These are the equations of sixteen coupled and forced har- 
monic oscillators, with forcing terms 5^^^ • (7^°) + C^°'> ■ S'^^^^ 
that are purely oscillating. Since C is real, the solutions are 
easily obtained in the form 

C(i)(t) = ci e*^"+c*i e"'"*. (A14) 



Recall that the matrix S{t) depends on e and 77 via the classical 
values and an additional explicit term —ujm ecos{ilt). If we 
make use of the approximate solution found before, we can 
thus identify a matrix S'^") independent of the perturbation, 
a matrix 5*^^^ linear in e and 77 and a matrix 5*^^^ quadratic 
in e and 77. Again, we look for a perturbative solution for the 
matrix C, i.e 



Hence, to first order in e and 77, the effect of the modulation on 
the correlations is to add an oscillating term with frequency 
il and mean value 0. As in the case of an unidimensional 
resonator, the amplitude of this oscillation will be greater 
when the modulation frequency fl is chosen in resonance 
with the eigenfrequencies of the normal modes. Finally, the 
equations for second order are 



These are again the equations of sixteen coupled and forced 
harmonic oscillators, but this time the forcing terms S"*^^' • 
C(i) + • and ^(2) . c-Co) + ^(o) . ^(2)t ^ave also 

a constant part. The solutions are easily obtained in the form 



3. Asymptotic behavior in the single modulation regime 



a. Classical solution. 



C3 e 



ant 



-i20t 



(A16) 



Hence, as for the linear solutions, to second order in e and rj, 
the effect of the modulation on the correlations is to add a con- 
stant shift and an additional oscillating term with frequency 
257 and mean value 0. 



We can come back to the single modulation scenario by 
putting 7; = in the above analysis. By doing so, we can get 
approximate analytic expressions for the asymptotic mean 
values A{t), Q{t) and P{t). However the complete formulas 
are too long to be reported here and we must limit ourselves 
to a semi-numeric expression, where we substitute all values 
as in Sec. Ill except for the interesting parameter e. For 
example we report the expression of the mirror position Q{t) 
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Q{t) = 14684.7 - 2784.43 

+ e(4947.11 cos{nt) - 14.79 sin(l]i)) 

+ (l64.97 cos(2m) - 0.50 sin(2m)) . (A17) 

As anticipated in the main text, to first order in e the mean 
values have an asymptotic oscillatory behavior, which well 
describes the exact asymptotic solution. To be precise 
however, we cannot neglect the second order contributions: 
indeed, while second harmonic oscillations are one order of 
magnitude smaller, the constant shift is comparable to first 
order effects and must be taken in account. 



We conclude this section with one last comment on why 
the number of phonons is constant in time. We know that 
the position Q{t) and momentum P{t) of the mirror oscillate 
with frequency ^ ojm and a relative phase shift of tt/2 (slight 
modifications being induced by the interaction with the optical 
subsystem). In the same way Cn = (Sq^) and C22 — (Sp^) 
oscillate with twice this frequency, i.e. ^ 2aj„j, and with twice 
this relative phase shift , i.e. tt. In turn, the two oscillations 
cancel each other out when summing Cn and C22, thus giving 
a time-independent number of phonons. In addition we see 
that the modulation is most effective on the mirror correlations 
when il ^ 2uj,n, as stated before. 

4. Asymptotic behavior in the two modulation regime 

a. Classical solution. 



b. Linearized quantum solution. 

We can also look at the quantum properties of the system 
in the asymptotic regime, as a function of the modulation 
strength e. Fixing all other parameters to values in the text, 
we find for example the following expression for the number 
of phonons in the mirror np;io„ (t) w (Cii(t) + C22(i) ~ 

nphon{t)^OM + e^ 4.14 

+ e(^0.14 cos(m)-0.01 sin(m)) 

- (0.02 cos{2flt) - 0.21 sin(2m)) . (A18) 

For completeness we also report the expression for the single 
correlation Cn 

Cii(0 = 0.56 + 6^ 4.01 

+ e(o.28 cos(m) - 1.63 sin(m)) 

+ (0.03 cos(2m) - 0.20 sin(2m)) . (A19) 

We see that Cii{t) (and similarly 6*22(0) has strong oscil- 
lations in time proportional to e. Hence we can say, at least 
qualitatively, that squeezing will be dominated by first order 
effects in the range of values considered. On the contrary 
the number of phonons can be considered time-independent, 
with oscillations that are negligible if compared to the con- 
stant term. Moreover, we know that the mirror is cooled close 
to its ground state when the system is unmodulated. There- 
fore the number of phonons is strongly dependent on e^, and 
the constant shift due to second order effects becomes quickly 
the dominant effect. These results agree very well with the 
numerical simulation summarized in Fig. [3] 



We now reintroduce the second modulation and study the 
interplay between the two. Again we would like to fix all 
parameters except e, 77 and (p to the values found in the main 
text. However, already at the classical level, expressions for 
A{t), Q{t) and P{t) tend to become rather long and complex 
since we have now 3 free parameters. Therefore we will 
substitute also the numerical values of e and 77 (values are 
found in Sec. |1V| |. This is not so bad: indeed recall that we 
are particularly interested in the dependence of quantum 
properties on the relative phase <j). For example we report the 
expression of the mirror position Q{t) 

Q{t) = 17523.4 - 357.13 co&{(j)) + 315.98 sin(0) 

+ 1484.13 cos(m) - 4.43 sin(rjt) 

+ 2201.22 cos(m + </>)- 1870.83 sin(m + cf) 

+ 14.84 cos(2m) - 0.04 sin(2rjt) 

+ 22.62 cos{2nt + 0) - 18.78 sin(2f7i + 0) 

+ 113.53 cos(2m + 20) - 32.43 sin(2m + 2(1)). 

(A20) 

Without losing much time on the cumbersome formula above, 
we only point out that again first order effects are dominating 
(as we can see by comparing oscillations at and oscillations 
at 2ri). However we also see that, already at the classical 
level, the the phase (f) has a strong influence on the amplitude 
of oscillations. This is a clear sign that the phase plays indeed 
an important role in the system dynamics. 



b. Linearized quantum solution 

We turn now to quantum properties of the system. Again 
we fix all parameters at the values used above, except for the 
relative phase 0. To second order in the perturbation, we get 
expressions for the number of phonons and the correlation 
Cii (t) that depends on the phase as 
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nphonit) = ^ (l.l67 + 0.087 cos(0) + 0.753 sm{(f>) + 0.086 cos(m) - 0.005 sin{nt) 

~ 0.004 cos(2rJi) + 0.037 sin(2m) + 0.008 cos(f7i + <p) ~ 0.006 sin(m + </)) 

- 0.024 cos(2r2t + + 0.005 sm(2r2t + 0) + 0(10"^)) , (A21) 

Cii(t) =1.09 + 0.02 cos(0) + 0.39 sin(0) + 0.08 cos(m) - 0.48 sm(m) 

+ 0.002 cos(2m) - 0.02 sin(2m) + 0.39 cos(m + 0.05 sm{nt + <j)) 

+ 0.001 cos(2Sl< + 0) - 0.002 sin(2m + (/>) + 0(10"^), (A22) 



where for brevity we have neglected the smallest terms. Look- 



ing at expression ( A21 1 above, it is clear that the main effects 
of the modulation are contained in the first and third terms, 
other terms being an order of magnitude smaller than the two. 
The two bigger terms are both independent of time (similarly 
to the single modulation case), hence they must come either 
from the unmodulated solution C^"' or from the constant part 
of the second order solution C^^\ Again, since in the unmod- 
ulated scenario the mirror is very close to its ground state, we 
can reasonably assume that the matrix C'"-* contributes in a 
negligible way. Therefore, the dependance of Umax on the 
phase (j) is almost entirely described by the matrix C'^-' and is 
a second order effect in the modulation strengths e and 77. The 
maximum number of phonons (i.e. Eq. ( A21[ l maximized over 
one period r — 27r/r2 of evolution), as well as the various con- 
tributions described here, are plotted in Fig. |6] We see that the 
analytic approximation correctly resembles the numeric solu- 
tion (see Fig. |4]i, the two differing only by a small constant 
shift which is due to higher order corrections. However, the 
qualitative behavior is fully understood already at the second 
order; therefore we do not report here explicitly higher orders 
contributions. 



From Eq. (A22i, we see that Cn (and similarly C22) un- 
dergoes strong oscillations in time, reaching a minimum value 
that depends much on the phase (j). This tells us that mechani- 
cal squeezing will have a very similar behavior and hence will 
also depend strongly on 



Entanglement and quantum discord have instead a much 
smaller response. Indeed, both quantities are computed using 
all entries of the matrix C (and not only two entries as in the 
case of phonons). Each entry will have an expression similar 
to (A22i and made of three parts: a constant part, a time 
independent part which oscillates with the relative phase (p 
and a time dependent part. This can be written in the general 
form 



a. 



A, 



(A23) 



with A. 
B 



cos( 



Bij and ipij constants. In general, the oscillations 
ij K.^oy:j) + (fiij) will be out of phase with one another. Also 
the time dependent parts (t) will be generally oscillating 
out of phase. Therefore, summing many entries together, 
these parts will cancel each other out (as a sum of incoherent 
waves) and only the constants Aij survive to play a relevant 
role. From this handwaving reason, we expect that entan- 
glement and quantum discord should be quite insensitive to 
phase ip, as is the case in Fig. |4] 




1 2 3 4 5 6 







Figure 6. In blue (upper solid) the maximum number of phonons, 
from Eq. ( |A21^ , plotted against the relative phase (f). Contributions 
due to different expansion orders are explicitly included: in black 
(lower solid) the number of phonons in the unmodulated case (from 
C'"'); in red (upper dashed) the second order time independent con- 
tribution (from C'^'); in green (lower dashed) the maximum over 
one period of the time dependent contributions (from first order C''^' 
and second order C'^'). The total (blue) curve is equal to the sum of 
the other three curves. 



